##make a biplot
source('/home/jc148322/scripts/libraries/cool_functions.r')
library(SDMTools)


pre.file='/home/jc148322/frog_declines/Pre_decline.csv'
post.file='/home/jc148322/frog_declines/Post_decline.csv'
climate.dir='/home/jc165798/Climate/CIAS/AWT/250m/baseline.76to05/bioclim/'; setwd(climate.dir)

bc01=read.asc.gz('bioclim_01.asc.gz')
bc12=read.asc.gz('bioclim_12.asc.gz')

pos=make.pos(bc01)

pre.occur=read.csv(pre.file); colnames(pre.occur)=c('spp','lat','lon')
post.occur=read.csv(post.file); colnames(post.occur)=c('spp','lat','lon')

vois=c('pos','pre.occur','post.occur')

for (voi in vois) {
	tdata=get(voi)
	tdata[,'bc01']=extract.data(cbind(tdata[,'lon'],tdata[,'lat']),bc01)
	tdata[,'bc12']=extract.data(cbind(tdata[,'lon'],tdata[,'lat']),bc12)
	if (voi==vois[1]) tdata=tdata[,c('bc01','bc12')] else tdata=tdata[,c('spp','bc01','bc12')]
	assign(voi,tdata)
}

#edit and simplify the data for biplot
pos[,1] = round(pos[,1],1);
pos[,2] = round(pos[,2]/10)*10 #round bc01 and bc12
pos = unique(pos) #find only unique combinations of rounded bc01 and bc12



out.dir='/home/jc148322/frog_declines/SDM/regular/spp_biplots/';dir.create(out.dir); setwd(out.dir)
#CREATE PLOTS AND IMAGES

species=c('dayi','nanno','serrata','rheoco')

for (spp in species) {
png(paste(spp,'.png',sep=''), width=8,height=8,units='cm', res=300,pointsize=5)

        plot(pos, xlab="Mean Annual Temperature",ylab='Annual rainfall (mL)',xlim=range(pos[,1],na.rm=T),ylim=range(pos[,2],na.rm=T), type='n',cex.lab=1, cex.axis=1, font.lab='2')
        points(pos, col='grey', pch=19,cex=3)
		
		predata=pre.occur[which(pre.occur$spp==spp),c(2:3)]
		postdata=post.occur[which(post.occur$spp==spp),c(2:3)]
		points(predata,col='black',bg='forestgreen',lwd=0.5,pch=21)
		points(postdata,col='black',bg='red',lwd=0.5,pch=21)
		
		# legend(23,7000, c('Climate space', 'New Bd SDM', 'Original Bd SDM'), fill=plotcols, cex=1, bty='n')
		 legend(23,7000, c(paste(spp,' - pre',sep=''), paste(spp,' - post',sep='')), col=c('black','black'),pt.bg=c('forestgreen','red'),pch=c(21,21), pt.lwd=0.5, cex=1, bty='n')
		
dev.off()
}
